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An understanding of how an initially Gaussian error volume becomes non-Gaussian over 
time is an important consideration for space-vehicle conjunction assessment. Traditional 
assumptions applied to the error volume artificially suppress the true non-Gaussian nature 
of the space-vehicle position uncertainties. For typical conjunction assessment objects, 
representation of the error volume by a state error covariance matrix in a Cartesian 
reference frame is a more significant limitation than is the assumption of linearized 
dynamics for propagating the error volume. In this study, the impact of each assumption is 
examined and isolated for each point in the volume. Limitations arising from representing 
the error volume in a Cartesian reference frame is corrected by employing a Monte Carlo 
approach to probability of collision (P c ), using equinoctial samples from the Cartesian 
position covariance at the time of closest approach (TCA) between the pair of space objects. 
A set of actual, higher risk (P c > 10 -4 ) conjunction events in various low-Earth orbits using 
Monte Carlo methods are analyzed. The impact of non-Gaussian error volumes on P c for 
these cases is minimal, even when the deviation from a Gaussian distribution is significant. 
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NTW 

xyz 

Vectors 


= right ascension of the ascending node 
= argument of perigee 
= true anomaly 

= Normal, in-track, cross-track Cartesian coordinates centered on space-vehicle 
= Encounter frame coordinates ( x in direction of relative position, y in direction of relative velocity) 

are indicated by bold, lower-case symbols. Matrices are indicated by bold, upper-case symbols. x = -^(x); 


x = unit vector of x ; \ T = transpose of x ; is the expectation operator; I is the identity matrix. 


I. Introduction 

C ONJUNCTION assessment risk analysis (CARA) is the process of monitoring and evulating the collision risk 
from unique close approach events between two space objects. The “big sky” theory, in which it is presumed 
that object densities are small enough that the likelihood of collisions between spacecraft is discountable, is no 
longer valid as evidenced by recent collisions . 1 Satellite protection is critical both to maintain the satellite mission 
and to avoid polluting the orbital neighborhood. 

Collision risk can be expressed in several ways, most commonly by the miss distance (the closest in position the 
two objects will come) or, better, by the probability of collision (P c ). The P c metric is a useful mechanism to 
encapsulate risk but it requires accurate knowledge of the state uncertainty, typically represented by the state error 
covariance matrix. The analytic calculations of P c rely on an assumption that the propagated state error distribution 
is Gaussian . 2 This Gaussian assumption appears to be on questionable ground, and the impact of this non-Gaussian 
behavior on P c is not well understood . 3 The orbit determination process results in an initially Gaussian error 
distribution, generally assuming that the observation residuals are independent and identically distributed from a 
zero-mean Gaussian and the dynamics are suitably modeled. In practice, however, these assumptions may not be 
entirely true. Assuming an initially Gaussian error volume, longer-term non-Gaussian behavior in the error volume 
often becomes evident after a couple of days of propagation, even for well-tracked objects . 4 Operations for the 
National Aeronautics and Space Administration (NASA) Robotic CARA task regularly require predictions a week 
or more in the future to determine potential conjunction assessment events of interest. Therefore, a significant 
deviation from a Gaussian error volume may be expected. Prior to developing a solution to mitigate or fix non- 
Gaussian effects, this investigation addresses the nature of non-Gaussian effects from propagation for low-Earth 
orbit (LEO) objects and, subsequently, determines the resulting impact to P c using a Monte Carlo approach on a set 
of actual, high risk events from the NASA Robotic CARA task. 

II. Background on Non-Gaussian Error Volumes 

An initially Gaussian error volume representing space-vehicle uncertainties will become non-Gaussian over 
time, which is a consequence of the nonlinear dynamics governing space-vehicle motion. This nonlinearity is a 
fundamental aspect of the dynamical system and its coordinate system representation . 5 To correctly process how an 
initial Gaussian error volume evolves over time, nonlinear dynamics must be utilized. The effects from the nonlinear 
dynamics on the error volumes becomes more evident for larger initial error volumes and for longer propagation 
times. 

This section first reviews how the true non-Gaussian nature of the error volume is artificially suppressed by 
using a linearized dynamical system to propagate the error volume (linearized dynamics) and by representing the 
error volume in a Cartesian reference frame. Then, there is a discussion of methods to recover, in part or in whole, 
the true non-Gaussian nature of the error volume. Finally, there is a review of research that relates non-Gaussian 
error volumes to conjunction assessment. 

A. Source of Evolution to Non-Gaussian Error Volumes 

The state error covariance matrix P is routinely used to represent the uncertainty of a satellite state vector x in 
terms of a Gaussian ellipsoid. This covariance may not match the true state uncertainty; if they do agree, Horwood et 
al . 6 refer to this condition as uncertainty consistency. In Junkins et al ., 7 an investigation of how an initially Gaussian 
error volume becomes non-Gaussian over time is presented. In that research, they present two aspects where 
linearization assumptions with regard to how the error volume is handled can suppress the true non-Gaussian nature 
of the error volume. The first aspect relates to the assumption of linearized dynamics in error theory, as the entire 
error volume is propagated based on a linearized Jacobian evaluated at the nominal state (refer to the Appendix for 
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additional information). The linearized dynamics assumption is summarized by the similarity transformation for 
propagating covariance from one time t 0 to another time t using the state error transition matrix, ®, that is, 


p(0 = ®(*> ? 0 )P(? 0 )® r (Mo)- (1) 

The second aspect where a linearization of the error volume appears is with transformations between reference 
frames. If one frame x relates to another frame y by 


y = g(x) , 


( 2 ) 


then the covariance in frame x can be transformed to frame y using a similarity transformation, that is, 

P y =AP x A r , 


(3) 


where 


A = 



(4) 


is the linearized Jacobian evaluated at the nominal state. 

A Cartesian representation for the covariance, as is used for current CARA operations, is an unfortunate choice 
for representing the error volume of a nonlinear trajectory. Similarity transformations for covariance, by their nature, 
are shape-preserving, so an initially Gaussian error volume will artificially remain Gaussian for all time and in any 
reference frame. If the error volume becomes too large, whether due to a long propagation time, a large initial 
uncertainty volume, or some combination of the two, then these linearization assumptions can result in the 
covariance being an inadequate representation of the actual uncertainty volume - a violation of uncertainty 
consistency. The non-Gaussian effects can become evident on a timescale of hours to days. This situation is most 
challenging when dealing with uncorrelated target (UCT) initial orbit determination (IOD) objects, as the initial 
uncertainty volumes are large and problematically- shaped, but non-Gaussian error volumes must be considered for 
conjunction assessment activities as well. 

B. Techniques to Recover the Non-Gaussian Error Volume 

Methods for restoring the true non-Gaussian state error distribution can be classified into several categories. The 
first is the familiar Monte Carlo strategy. The second involves the selection of an alternative coordinate 
representation for the error volume to address the limitations of the Cartesian reference frame, and the third 
addresses the linearized dynamics assumption for propagating the error volume. 

The Monte Carlo method, taking pseudo-random Gaussian perturbations of the state uncertainty at epoch and 
propagating these perturbed states, is an accurate representation of the error volume and is often employed as a truth 
source for the propagated uncertainty state. In representing the error uncertainty as a cloud of points that are 
propagated with the full nonlinear dynamics, there are no linearization assumptions involved. Monte Carlo methods 
can come with a severe computational cost as all perturbed states must be propagated, which has been a motivational 
factor to formulate alternative approaches. 

Another technique to recover the non-Gaussian error volume is to select an alternative set of coordinates in 
which the covariance more closely approximates the true error volume. One approach is to change the coordinates 
used for propagating the covariance, such as classical orbit elements, 8 equinoctial orbit elements, 4,7 ’ 910 Poincare orbit 
elements, 11 ’ 12 polar coordinates, 7 and curvilinear coordinates. 4 ’ 13 ’ 14 Another approach is to simply transform the 
Cartesian covariance into a more suitable representation using Eq. (3) at the end of the propagation time. Sabol et 
al. 4 demonstrate this method by transforming the covariance from osculating Cartesian Earth-centered inertial (ECI) 
into osculating equinoctial. Both approaches to this technique are computationally efficient, as only a single state 
vector and the associated covariance matrix are propagated. 

Addressing the linearized dynamics assumption for propagating the error volume is challenging as this 
assumption relates to how the state uncertainty is propagated. One approach is to replace the linearized dynamics for 
propagating the error volume with a higher order error theory, generalizing the state error transition matrix as a state 
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error transition tensor. 11,15,16 Another approach employing Gaussian mixture models 17 addresses the limitations of 
representing the error volume as a single Gaussian probability density function (pdf). Gaussian sum filters estimate 
any arbitrary pdf as a weighted sum of Gaussian components and have been successfully applied to determining 
satellite error volumes. 10,18,19 An added benefit is that Gaussian sum filters may make use of existing uncertainty 
mapping techniques, such as those employed by extended Kalman filters or unscented Kalman filters, for 
propagating the Gaussian error volume components. One can also represent the error volume pdf with other 
formulations, as demonstrated by Edgeworth filters 6 and polynomial chaos expansions. 12 This collection of methods 
shows potential and could be competitive or superior to a Monte Carlo approach. 

Techniques that alter the coordinate representation for the error volume and address the linearized dynamics 
assumption for propagating the error volume are at times intertwined. Junkins et al. 7 found that an equinoctial 
formulation and propagation of covariance performed better than polar coordinate formulation, which in turn 
performed better than Cartesian coordinate formulation, although the authors did not separate the two sources of 
nonlinearity. Sabol et al. 4 separate the impact of linearized dynamics from the Cartesian comparison frame in their 
seminal work. Hoots 8 maintains that an approximate formulation for the error volume, represented analytically as a 
polynomial expansion in classical elements, eliminates nonlinearities in both the coordinates and in the propagation 
of the error volume. Folcik et al. 9 claim that mean equinoctial propagation of covariance performs better than 
osculating equinoctial propagation in containing the truth state. Even a subtle change, replacing the semi-major axis 
by the mean motion in equinoctial orbit elements, can result in better agreement with the true error volume. 10 

C. Non-Gaussian Effects and Conjunction Assessment 

While the theory of mitigating non-Gaussian error volumes is relatively mature, the impact of non-Gaussian 
error volumes to conjunction assessment is not well-characterized at present. Tanygin and Coppola 20 assess the 
impact of the bias term from a second-order error volume propagation on an entire space object catalogue with 
synthetic covariance values. They examine differences in time of closest approach (TCA), miss distance, and 
maximum P c associated with conjunction events and conclude that this second order bias effect can be significant, 
e.g., for an object about to decay. The second order bias addresses only one aspect of the impact of non-Gaussian 
error volumes, as an underlying Gaussian assumption remains and no attempt is made to address the limitations of 
representing the error volume in Cartesian coordinates. 

In their work, Sabol et al. 21 examine the impact of non-Gaussian error volumes on the miss distance distribution. 
In their Monte Carlo framework, they perform perturbations of the initial covariance matrix and compare the results 
with both the standard analytic formulation and with Monte Carlo perturbations taken at TCA. Their analysis reveals 
differences in the distribution of miss distances that are attributable to non-Gaussian error volumes. However, they 
make no attempt to demonstrate resulting changes in P c . The P c metric is the integral of the miss distance pdf with 
limits of integration from zero up to the combined physical size of the objects. Therefore, the only changes in the 
miss distance distribution that are significant to P c are those corresponding to small miss distances. Understanding 
the impact to P c from non-Gaussian error volumes is the most fundamental aspect from a CARA perspective and is 
the objective of this study. 


III. Propagated Error Volume Investigations 

Prior to investigating the impact of a non-Gaussian error volume on conjunction assessment events, this section 
explores the propagation of an initially Gaussian error volume by considering both propagated Monte Carlo samples 
and the propagated covariance. Then follows a discussion of representing the propagated error volume as equinoctial 
samples of the propagated Cartesian covariance, which is a mechanism to address the limitations of the Cartesian 
reference frame. Finally, the errors associated with the assumption of linearized dynamics for propagating the error 
volume are isolated and found to be negligible for typical conjunction assessment objects. 

A. Position Uncertainty Misalignment Angle 

For the initial investigation into propagated error volumes, the time evolution of an initially Gaussian error 
volume for a single object in a FEO orbit is considered for objects ranging from 400 km to 1300 km in altitude. A 
high-fidelity model in FreeFlyer® software is used in the propagation of the nominal and perturbed state vectors; the 
nominal covariance is propagated numerically using the FreeFlyer® state error transition matrix function. These 
Monte Carlo perturbations are generated at the initial epoch and represent the true error volume when propagated 
forward in time. 

It is instructive to go through one example in detail as its results are representative of the other scenarios. The 
initial Keplerian state { a , e , /, 12, m \ v} is {7718.4 km, 0.00039, 66.0°, 282.7°, 100.5°, 52.1°} and the initial ECI 
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position covariance diagonal elements are 15.7 m, 3.9 m, and 6.6 m respectively. Figure 1 displays the results from 
propagating the perturbed states for ten days in addition to the propagated Cartesian covariance. The coordinates are 
relative to the nominal state, located at (0,0) in the plot. The convention for representing the satellite state vector, 

x = [r vf where r is the position vector and v is the velocity vector, is the NTW Cartesian reference frame 
which is defined as 


* v * r x v * 

T = -j— r; W = t N = Tx W . (5) 

v rxv 

The AT-plane represented in Fig. 1 depicts the nominal orbit plane. The propagated Monte Carlo samples that 
represent the error volume deviate from a Gaussian distribution. Note that the motion of the space-vehicle is 
necessarily in the T -direction. The distribution of the error volume in the out-of-plane IT-component appears to be 
Gaussian. 



Figure 1. Monte Carlo perturbations from nominal state, propagated for 10 days, with the corresponding 3 o 
Cartesian covariance ellipsoid and the position uncertainty misalignment angle. The axes are not on a 1:1 
scale to emphasize the non-Gaussian distribution of the propagated Monte Carlo states. 

From Fig. 1 it is clear that the major principal axis of the position covariance (i.e., the red ellipsoid) is not 
perfectly aligned with the nominal velocity vector, an attribute that has been noted previously by Vallado and 
Alfano. 14 Here we define the position uncertainty misalignment (PUM) angle to be the angle between the major 
principal axis of the position covariance in the AT-plane and the T-axis as determined by the nominal state. This 
PUM angle varies over time, as shown in Fig. 2, dampening in an asymptotic manner over long propagation times. 
This position misalignment angle and its time variability is an attribute associated with propagation of covariance 
from using a state error transition matrix calculated from the osculating dynamical equations of motion. 

Figure 3 displays the same propagated Monte Carlo samples as shown in Fig.l with the addition of the nominal 
orbit trajectory. The propagated error volume is not aligned with the nominal orbit trajectory. Using the PUM angle, 
one could rotate all of the Monte Carlo values by this angle, or, equivalently, rotate the orbit trajectory in the 
opposite way, as represented in Fig. 3. Performing this rotation results in visual agreement between the propagated 
Monte Carlo samples and the rotated orbit trajectory. Although not shown here, an analogous treatment can be 
performed for the velocity components, where the velocity misalignment angle represents the difference between the 
major principal axis of the velocity covariance in the AT-plane and the nominal acceleration vector. 

Although Fig. 3 is simply a visual representation, this suggests that representing an osculating ECI covariance as 
an osculating orbital-based covariance may be successful in recovering the fundamental non-Gaussian error volume. 
This technique will be pursued in the next subsection. As an aside, the contribution from the osculating covariance 
propagation (as evident by the position misalignment angle) is of the same order of magnitude as the non-Gaussian 
behavior for this example. Some researchers advocate mean equinoctial propagation 9,22 of covariance to reduce the 
linearized dynamics errors, but by doing so there is the possibility of losing out on a genuine characteristic of 
propagation that could be important, especially for shorter propagation times. 
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Figure 2. Position uncertainty misalignment (PUM) 
angle (log scale) in terms of propagation time. 



’° %0 AO -30 -20 -10 0 10 20 30 40 50 

Relative Position - Tangential Component, km 

Figure 3. Same Monte Carlo samples as Fig. 1 with 
addition of the nominal and rotated orbit trajectory. 
The PUM angle is used for the rotation. 


B. Sufficiency of Addressing Cartesian Coordinate Limitations 

This technique to address Cartesian coordinate limitations is the Jacobian methodology proposed by Sabol et al. 4 
The osculating ECI Cartesian covariance at the end of the propagation time is converted using a Jacobian matrix into 
an osculating equinoctial covariance. Gaussian samples are computed from this osculating equinoctial covariance. 
Rather than remaining in equinoctial element space, these samples are converted back to Cartesian ECI using the 
full nonlinear transformation so that they can be compared directly with the propagated Monte Carlo samples. 
Figure 4 shows the comparison in the nominal orbit plane for both position and velocity respectively. 

The visual agreement in Fig. 4 is striking, indicating that addressing the Cartesian frame limitations is adequate 
for this example. This scenario is not unusual, and it is consistent with findings by Sabol et al. 4 and Horwood et al. 10 
For Sabol’ s “catalog-class scenario,” converting a covariance from osculating Cartesian into osculating equinoctial 
restores a distribution that is consistent with a Gaussian distribution even after ten days of propagation. Horwood 
indicates that accurate techniques such as Gaussian sum filters may not be necessary in some situations (e.g., small 
initial uncertainty in semi-major axis). Their emphasis is on the UCT IOD problem and the smallest initial 
uncertainty in semi-major axis considered for their study is 1 km, which is several orders of magnitude greater than 
is typical for conjunction assessment situations. 




Relative Velocity - T angential Component, m/s 


a) 


b) 


Figure 4. Equinoctial samples of Cartesian covariance compared with propagated Monte Carlo samples in 
terms of a) position and b) velocity. 


6 

American Institute of Aeronautics and Astronautics 


C. Isolation of Linearized Dynamics Errors 

As the technique in the previous subsection only addresses the limitations of the Cartesian coordinate 
representation, the purpose of this subsection is to attempt to isolate the other error source: the linearized dynamics 
errors. Instead of comparing entire distributions (such as by the expected Mahalanobis distance distribution 4 or by 
statistical goodness-of-fit tests), an alternative approach is to compare single points of the distribution and identify 
the linearized dynamics errors in each pair. 

This process for isolating errors is outlined in Fig. 5 using x and x as the nominal and perturbed state vectors. 
This process can be repeated for additional Monte Carlo perturbed states to generate a profile for the entire 
distribution. Consider a perturbed state vector x at the initial epoch t 0 

x (h))= x fe) + ^ x > (6) 

where Sx is a Monte Carlo Gaussian draw from the initial covariance. The perturbed state x(/ 0 ) , the nominal state 
x(t 0 ) , and the covariance corresponding to the nominal state are all propagated to a later epoch t. The state error 
transition matrix <d(^, t 0 ) is calculated and retained so that the position difference assuming linearized dynamics can 
be recovered. 

x(r) 



Figure 5. Methodology to isolate linearized dynamics errors from errors associated with Cartesian frame 
limitations. 

The position difference at the later epoch under the assumption of linearized dynamics is given by 

* LinDyn (0 = *0 X x (*0 ) “ X ( ? 0 )) • (?) 

As both initial state vectors are propagated, the actual state deviations are calculated: 

^actual W= X ( ? )- X (T (8) 

The difference between Eq. (7) and Eq. (8) is the deviation from a Gaussian distribution. Part of this deviation is 
attributable to Cartesian coordinate limitations and the remainder is attributable to the linearized dynamics. Just as 
the osculating Cartesian covariance is converted into osculating equinoctial covariance in the previous subsection to 
address the Cartesian coordinate limitations, the same technique can be applied to individual state deviations. Then 

determine the state deviation £x* that would result in the actual state deviation £x actual after application of the 
Jacobian A, as represented by 


<%* — A -> &. actual . ( 9 ) 

The difference between these two Sx’s in Eq. (9) is the error associated with Cartesian frame limitations. The 
remainder of the deviation is attributed to the linearized dynamics error and is calculated as 
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^‘-^LinDyn. ( 10 ) 

Two examples are considered. Both are synthetic examples with larger-than-typical covariance values in an 
attempt to try to highlight the linearization errors. The first case is a 450 km by 6200 km orbit with initial ECI 
position covariance diagonal elements of 0.2 km, 0.2 km, 0.32 km respectively, and a propagation time of two days. 
The non-Gaussian root-sum-square (RSS) position error associated with Cartesian coordinates and linearized 
dynamics appears in Fig. 6. The total non-Gaussian deviation is the sum of these two terms. The color bar represents 
the position error on a logarithmic scale. Each point represents the comparison of a single Monte Carlo perturbation 
as compared with the nominal state. The AT-coordinates have been modified in two ways to make the figure clearer. 
First, the nominal AT-coordinates are rotated by the position uncertainty misalignment angle so that the Cartesian 
coordinates error is aligned with the rotated T-axis. Second, as the non-Gaussian error magnitude is already 
represented by the color scale, the deviation from a Gaussian distribution is removed from each point in the 
displayed AT-coordinates. Thus the distribution of points is consistent with a Gaussian. For this example, essentially 
all of the linearized dynamics errors are in the tangential direction, which makes it more difficult to observe as the 
position covariance is largest in this direction. The second case, appearing in Fig. 7, is based on a 450 km circular 
orbit with initial ECI position covariance diagonal elements of 0.1 km, 0.1 km, 0.2 km respectively, and a 
propagation time of 12 days. Here the linearized dynamics errors have non-zero components in both the tangential 
and normal directions. It is interesting to note that Fujimoto et al. 11 observe that the linearized dynamics errors result 
in a bias in longitude, with the largest contribution to the bias coming from the ends of the error volume distribution. 
For a smaller-sized covariance, such as in Fig. 6, this bias in longitude would principally be observed in the 
tangential component, whereas in a larger covariance that covers a significant portion of the orbit, this bias in 
longitude would have both normal and tangential components. 



Figure 6. Color represents non-Gaussian RSS position errors from 2-day propagation of 450 km by 6200 km 
orbit due to a) Cartesian frame limitations, and b) linearized dynamics. 

In the examples depicted in Fig. 6 and Fig. 7, the linearized dynamics errors are small compared with the 
Cartesian frame limitations. Any numerical errors, such as those arising from the state error transition matrix and the 
Jacobian transformation, are included in this linearized dynamics error term. It is improbable that the numerical 
error vector would negate the linearized dynamics error vector for every point evaluated, so the general behavior is 
indicative of the linearized dynamics errors. Using a realistic (i.e., small) initial covariance results in negligible 
linearized dynamics errors. This study uses an osculating ECI propagation of the covariance whereas, operationally, 
the NASA Robotic CARA task exploits the results from Air Force Astrodynamics Support Workstation (ASW) 
propagation, which itself employs an osculating equinoctial propagation of the covariance. Equinoctial propagation 
of covariance has been shown to result in smaller linearized dynamics errors than Cartesian propagation of 
covariance. 7,10 Therefore, these results should be considered to be an upper bound on the linearized dynamics errors 
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Figure 7. Color represents non-Gaussian RSS position errors from 12-day propagation of 450 km circular 
orbit due to a) Cartesian frame limitations, and b) linearized dynamics. 

observed operationally. Based on their relatively small size, the effects from linearizing the dynamics for covariance 
propagation are ignored in the next section. 


IV. Impact of Non-Gaussian Error Volumes on P c 

This section investigates the impact of non-Gaussian error volumes on the P c calculation for conjunction 
assessment. First the Monte Carlo technique for this P c comparison will be presented, followed by the results by 
applying this technique to actual conjunction assessment events. 

A. P c Comparison Technique 

The Monte Carlo approach for calculating P c in a conjunction assessment involves perturbing both satellites in 
the conjunction event from their nominal state using pseudo-random Gaussian draws in a manner consistent with 
their calculated uncertainty volumes and then propagating the perturbed states to TCA. 23 If the miss distance at TCA 
is within the combined spatial extent of the satellites (typically expressed as a hard body radius), then this is counted 
as a collision “hit.” With x hits observed in n Monte Carlo trials, the point estimate for P c is simply 


n 

Perturbations are performed either at an initial epoch or at the TCA associated with the nominal spacecraft 
trajectories. Taking Monte Carlo perturbations at the initial epoch is the more accurate method as there is no 
assumption of linearized dynamics for propagating the error volume, but this practice is computationally intensive 
because of the long propagation times. Perturbing at TCA is a faster calculation, as the propagation time between the 
original TCA and the updated TCA is minimal. Furthermore, high fidelity propagation is not necessary when 
perturbing at TCA; at this scale, propagations employing only two-body dynamics agree with the full special 
perturbations (SP)-based approach for short propagation times. 21 Perturbing at TCA has the disadvantage in that 
linearized dynamics errors are inherent in the error volume, but these linearized dynamics errors are generally small 
and considered negligible. 

This investigation uses Air Force-provided initial states and covariances, as captured in the vector/covariance 
messages (VCMs), with subsequent propagation to TCA using SP astro standards software. 24 A total of 10 7 
perturbations are made for each satellite at TCA using two different methods: 1) Gaussian samples of the 3^3 
Cartesian ECI position covariance; and, 2) Gaussian samples of the 6x6 equinoctial covariance. Samples of the 6x6 
Cartesian ECI covariance are not used, because the condition number for the covariance matrix approaches the 
accuracy limit for double precision numerical calculations in some cases. For the short propagation times involved 
from the original TCA, the effects of the velocity terms in the 6x6 Cartesian ECI covariance are very small, 
indicating that little additional error is introduced by not considering them. Subsequent propagation to the updated 
TCA is performed using two low-fidelity techniques: linear relative motion and simplified dynamics (4x4 
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geopotential with drag and solar radiation pressure disabled). A combined hard body radius of 20 m is employed for 
determining whether or not the Monte Carlo perturbations result in a hit. 

The input data set for this study comprise 248 historical, high risk (P c > 10” 4 ) conjunction events in various LEO 
orbits from the NASA Robotic CARA task database, ranging from 2006 to the present. Orbital properties for the 
primary (asset) and secondary objects as well as the corresponding distribution of propagation times appear in Fig. 
8. Propagation times for the primary objects encompass the typical operational timeline for these conjunction 
assessments, or up to seven days. The secondary objects have a wider distribution of propagation times, because 
many of these objects are not as well-tracked as the primary objects. 



Perigee Height, km Propagation Time, Days 

a) b) 

Figure 8. Characteristics of primary and secondary objects used in P c comparison study, a) Distribution of 
perigee and apogee heights, and b) distribution of propagation times. 


B. P c Comparison Results 

For each conjunction event, the two Monte Carlo P c point estimates are calculated using Eq. (11) from sampling 
the 3x3 Cartesian ECI position covariance (“Cart”) and the 6x6 equinoctial covariance (“eq”) respectively, that is, 


PCart ~ 



X e_q_ 
n 


( 12 ) 


The combined results from the P c comparison are displayed in Fig. 9, which displays yp eq — p Cart )/ 'Peart > 

percent difference between the two Monte Carlo calculations. Each point in the plot corresponding to a single event. 
The agreement between the two Monte Carlo methods is within ±15% for all cases. From the perspective of the 
CARA task, an order of magnitude or more change in P c would be considered significant, so none of the differences 
observed here are of operational significance. 

Because Monte Carlo methods are based on random sampling, it is natural to question how much of the 
variability in Fig. 9 may be explained by the Monte Carlo samples. To answer this question, a statistical test is 
employed to compare two proportions; this test is separately performed for each conjunction event to compare the 
pair of point estimates calculated using Eq. (12). The null hypothesis for each statistical test is that the two P c values 
are equivalent: 


#0 : Peq ~ PCart ■ 


( 13 ) 
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Figure 9. Percent difference in Monte Carlo P c values between equinoctial samples and Cartesian samples. 
Each point represents one conjunction event. 

The samples must be independent for this hypothesis test. Independence is guaranteed by using each Monte Carlo 
perturbed state only once. Some other researchers “recycle” Monte Carlo samples for computational efficiency. 21,25 
However, samples were not reused here in order to remove any possible violation of the independence assumption. 

The number of hits in a proportion is governed by a binomial distribution, and Fisher’s exact test statistic 26 uses 
a hypergeometric function to evaluate the validity of the null hypothesis. If the p-vahie associated with Fisher’s test 
is either too large or too small, then the null hypothesis is rejected and there is consequentially a statistically 
significant difference between the equinoctial and Cartesian P c values. A % 2 test statistic or the Score test statistic 
using a pooled estimate may be employed instead of Fisher’s test, but these methods assume a normal 
approximation to the binomial distribution, and the hypergeometric calculation for this study is not unwieldy here 
given the modest number of hits. 

Some number of failures is to be expected when running Fisher’s test 248 times even if the null hypothesis were 
always true. Only six cases (out of 248) failed at the 95% level and two cases failed at the 99% level. The details of 
these failures appear in Table 1. Failure at the 95% level means that the />-value is either greater than 0.975 or less 
than 0.025. All failed cases were repeated using 4 x 10 7 trials, and only the third case remained as a consistent and 
extreme outlier. Here, “extreme” is in reference to the /?-value being small (3.0 x 10” 9 ); the two resulting P c values 
are still very similar: p eq = 2.77 x 10” 4 and Pc ar t = 2-56 x 10” 4 . This indicates that, with one possible exception, all 

of the variation in Fig. 9 is explainable by the Monte Carlo sample variations. The “trumpet”- shaped nature of the 
scatter plot is because of less statistical confidence at lower P c values from fewer Monte Carlo hits. Differences 
between the two methods could be exposed if a higher number of trials are run, but any actual differences in the P c 
values is necessarily small. 

T a ble 1. Events that originally failed at 95% level and subsequent results from repeated t est. 


Failed 
Event # 

Original Test ( n = 10 7 ) 

Repeated Test (n = 4 x 10 7 ) 

%eq 

%Cart 

Fisher p - value 

%eq 

%Cart 

Fisher p - value 

1 

1,154 

1,282 

0.996 

4,847 

4,938 

0.824 

2 

2,600 

2,449 

0.017 

10,219 

10,169 

0.366 

3 

2,800 

2,522 

7.3 x 10“ 5 

11,099 

10,248 

3.0 x 10“ 9 

4 

34,997 

34,427 

0.015 

138,816 

139,220 

0.779 

5 

50,820 

51,526 

0.987 

206,362 

205,052 

0.020 

6 

91,140 

92,048 

0.984 

366,963 

365,417 

0.035 
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The outcome of the P c comparison test is surprising, as it was expected going into this study that there might be 
a significant difference in the calculated P c value resulting from non-Gaussian error volumes. One explanation for 
these observed results would be that the actual deviations from a Gaussian error volume are insignificant. To 
estimate the non-Gaussian deviation for each object in each event, a point is selected corresponding to the 30- 
uncertainty value along the largest principal axis of the position Cartesian covariance matrix. Using the Jacobian 
technique described earlier, the non-Gaussian deviation due to Cartesian coordinates is calculated by comparing the 
original Cartesian position with the linearized position of the corresponding equinoctial point. The results appear in 
Fig. 10 for all 248 cases (primary and secondary objects) using a log-log scale. Many of the events correspond with 
a small deviation from a Gaussian distribution. The median non-Gaussian deviation (at 3cr) for secondary objects is 
on the order of meters, although the deviation exceeds kilometers in some cases. It is not surprising that small non- 
Gaussian deviations result in a negligible impact to P c , but it is unexpected that this also holds true for large non- 
Gaussian deviations. The single case that failed has a modest non-Gaussian deviation of 28 m for the secondary 
object. The data in Fig. 10 are consistent with a quadratic relationship (corresponding to a slope of two in a log-log 
plot) that can be understood analytically. This relationship means that the deviation of an error volume from a 
Gaussian distribution due to limitations of the Cartesian reference frame can be estimated from the largest sigma 
value of the diagonalized position covariance. 



Largest Sigma Value of Diagonalized Position Covariance, km 


Figure 10. Distribution of non-Gaussian position errors for conjunction assessment events in P c comparison 
study. 

Conjunction events can be visualized at TCA by using the encounter frame. 2 In this representation, x is in the 
direction of the relative position vector (between the primary and secondary objects) and y is in the direction of the 
relative velocity vector. Two conjunction events, the extreme outlier event and a conjunction event displaying a 
large non-Gaussian deviation for the secondary object, are represented in Fig. 11. For each event, the secondary 
object’s nominal state is located at the origin and the primary object’s nominal state is located on the x-axis at the 
miss distance. The state uncertainty for both objects is represented in terms of the 3 a Cartesian covariance ellipsoid 
as well as Monte Carlo equinoctial samples of the Cartesian covariance. In Fig. 11a, there is a slight deviation from 
a Gaussian distribution bending towards the negative z-axis for the secondary object which is consistent with a 
slightly higher value of the equinoctial P c . In Fig 1 lb, the significant deviation from a Gaussian distribution for the 
secondary object results in a statistically insignificant difference of only 1.8% in the Monte Carlo P c values. This 
case is interesting for several reasons. Most of secondary object’s position uncertainty is in the y -direction and, 
therefore, is not apparent in Fig. lib. The encounter duration for this event, as represented by the distribution of 


*Using trigonometric relationships, if the covariance is much smaller than the Earth radius, then the leading term for 
the non-Gaussian deviation from an in-track uncertainty in a circular orbit at 3 a is (3cr) 2 /la . In the sample set 
analyzed in this research, both the variation in a and deviation from e = 0 are small, resulting in the good agreement. 
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Figure 11. Encounter frame representation for two conjunction events: a) P c extreme outlier event, and b) 
an event corresponding to a large deviation from a Gaussian distribution for the secondary object. 

TCA times resulting from the Monte Carlo process, is approximately half of a minute. The long duration of this 
conjunction event and the non-Gaussian deviations are coupled, limiting the utility of the encounter frame for 
understanding this event. All of the examples with significant non-Gaussian deviations examined in this study result 
in P c values comparable to what would be expected for typical short duration encounters that satisfy Gaussian and 
rectilinear assumptions. The theoretical basis for this behavior is a topic for future investigations. 


V. Conjunction Assessment Operational Considerations 

The first step in conjunction assessment operations relies on the Joint Space Operations Center (JSpOC) at 
Vandenberg Air Force Base “screening” assets against other objects in the high accuracy space object catalog. This 
screening entails placing a keep-out volume (called a screening volume) around a primary object and predicting 
when another object will enter that volume. These volumes are determined by using an empirically based covariance 
approach 27 or by some other means, but are intended to be sufficiently large to capture the possible conjunctions of 
interest. This screening is based on a miss distance approach between the nominal states, and therefore this 
calculated miss distance is not affected by non-Gaussian deviations. Provided the screening volume is large in 
dimensions compared to the non-Gaussian deviations of the error volume, the potential conjunctions of interest will 
be captured. Then the true collision risk associated with these events is evaluated. All potential conjunction events 
are evaluated for risk, including an analytic P c calculation, and higher risk events warrant additional monitoring and 
activities. 

Based on the findings of this study for LEO events, if the P c is sufficiently high, then non-Gaussian deviations 
result in a negligible impact to P c and do not change the risk of the conjunction event. As lower P c events (P c < 10” 
4 ) are outside the scope of this investigation, a possibility exists that a conjunction event that is considered to have a 
low collision risk might have a significantly higher P c than expected because of non-Gaussian deviations in the error 
volume. This scenario can be mitigated by a simple P c Monte Carlo calculation that can be performed using Air 
Force conjunction assessment data products (whether the product is the Orbital Conjunction Message (OCM) 
currently available or the Conjunction Data Message (CDM) available in the near future) which capture the two 
objects’ state vectors and position covariances at TCA. Cartesian perturbations of the 3^3 position covariance can be 
transformed using the Jacobian, which accounts for the non-Gaussian deviations arising from the Cartesian 
coordinate limitations. A Monte Carlo approach to P c using low- fidelity propagation from the nominal TCA and a 
reasonable number (1 0 5 ) of Monte Carlo trials executes quickly and, although the accuracy of the Monte Carlo 
method is low, it is of sufficient quality to establish that an ultra-low P c value does not become unacceptably high 
because of non-Gaussian behavior. 


13 

American Institute of Aeronautics and Astronautics 


VI. Conclusions 

An initially Gaussian error volume will become non-Gaussian over time due to the nonlinear dynamics inherent 
in propagation. The combination of linearizing the dynamical system for propagating the error volume (linearized 
dynamics) and representing the error volume in a Cartesian frame artificially suppresses the true non-Gaussian 
nature of the error volume as expressed by the state error covariance matrix. The limitations from representing the 
error volume in Cartesian coordinates dominate over linearized dynamics errors for typical conjunction assessment 
objects. Neglecting the linearized dynamics errors simplifies the processing for conjunction assessment, as retention 
of existing theories and propagation methods is permitted. A sufficient method to account for the non-Gaussian 
behavior addresses the limitations of Cartesian coordinates in representing the error volume: convert Gaussian 
Monte Carlo samples from the Cartesian covariance at TCA using the linearized Jacobian into osculating equinoctial 
elements, and use these perturbed states as the basis for a Monte Carlo calculation of P c . Furthermore, performing 
this technique on a set of 248 LEO conjunctions with P c > 10” 4 results in an operationally insignificant impact to 
calculating P c , even in cases where the non-Gaussian deviation is significant. Although this is not an exhaustive 
sample, the results suggest that the impact to P c from non-Gaussian error volumes is insignificant for high P c cases. 
As state uncertainties represented by the covariance are generally small, and if the covariances are small, the non- 
Gaussian deviations are consequently small. The impact of non-Gaussian error volumes on low P c events is not 
examined in this study but is mitigable using data products currently provided by the JSpOC using the proposed P c 
Monte Carlo framework. 


Appendix 

Let the spacecraft state vector equations of motion be represented by 

x = f(x,^). (Al) 

State deviations at an initial time t 0 can be mapped to state deviations at another time t using a state error transition 
matrix <D: 


Sx(t)=®(t,t 0 )Sx(t 0 ). 

The state transition matrix can be numerically integrated by considering the following equation: 

<b(M 0 )=A(f)®(M 0 ) 


subject to the initial condition 


®(Wo) = I. 


where I is an identity matrix of order six and the Jacobian matrix A(/) is defined by 


dx(t) 

dx(t) 

dx(t) dx(t) 


dx(t) 

dy(t) 

dz(t ) dz(t) 


m) 

W) 

W) 

ll 

■+0 x- V 

56 ^ 

II 

dx(t) 

dy(t ) 

dz(t ) 

dz(t) 

dx(t) 

8z{t) 

dy(t) 

dz(t) . 

dz(t) 


dz(t ) 


dz{t) 


dx(t) 


dz(t ) 


(A2) 


(A3) 


(A4) 


(A5) 


This notation for this Jacobian presumes osculating ECI propagation as is used in this paper. The covariance matrix 
at time t can be expressed as 
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P(/) = J £ , |^x(/)^x r (/)| 


(A6) 


where is {•} is the expectation operator. Substituting Eq. (A2) into Eq. (A6) leads to the result that covariance 
propagation can be written as 


p (0 = °(MoM*o)^ r (Mo), (A7) 

which is the same similarity transformation shown as presented in Eq. (1). The coordinate representation of the 
dynamics is an inherent aspect of the nonlinearity of the problem, so changing the reference frame of propagation in 
Eq. (Al), for example from osculating ECI to osculating equinoctial, will affect the resulting linearized dynamics 
errors. 
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